function out = display_image_fun(datapath)

% y = h5read(datapath,'/Image/image');
% y = double(y); % /!\ int16 not automatically converted into double

[pX,pY,OD] = ProcessImage(datapath);

% Estimate number of atoms using camera
amp = (pX(1) + pY(1))/2; % Take mean amplitude of two Gaussian fits.
capture_fraction = 0.5*(1 - sqrt(1 - (expp.na_telec)^2)); % fraction of total light captured by lens ( ~ 0.5 %)
total_num_electrons = amp*sqrt(2*pi)*mean([pX(3) pY(3)]);
total_photons = total_num_electrons/(expp.qe_std*expp.trans780*expp.trans_telec*capture_fraction);
num_atoms = 2*total_photons/(img_time*cons.Gamma) % num_atoms = (num_photons/time) / (scattering rate)

figure(1)
subplot(3,3,[1 2 4 5])
imagesc(OD);
if ROI
    line([1 1392],[ROI(1)+wndw ROI(1)+wndw])
    line([1 1392],[ROI(1)-wndw ROI(1)-wndw])
    line([ROI(2)+wndw ROI(2)+wndw],[1 1040])
    line([ROI(2)-wndw ROI(2)-wndw],[1 1040])
end
axis equal
xlims = xlim;
ylims = ylim;

subplot(3,3,[3 6])
plot(yInt,yC,gauss(pY,yC),yC);
set(gca,'YDir','reverse');
ylim(ylims);

subplot(3,3,[7 8])
plot(xC,xInt,xC,gauss(pX,xC));
xlim(xlims);

% Give the width as the std deviation in um
px_size = 6.45; % um
mag = 1/2; % magnification
pY(3) = (pY(3))*px_size/mag;
pX(3) = (pX(3))*px_size/mag;

h = subplot(3,3,9,'Visible', 'off');
cla(h);
text(0, .5, sprintf('%s\n%s\n%s\n%s', ...
    ['v-width (um): ' num2str(pX(3))], ['h-width (um): ' num2str(pY(3))], ...
    ['peak(v,h) = (' num2str(round(pX(2))) ',' num2str(round(pY(2))) ')'], ...
    ['Na = ' num2str(num_atoms,'%e')]), ...
    'VerticalAlignment', 'middle', ...
    'FontName', 'Courier New', ...
    'FontWeight', 'bold', ...
    'FontSize', 20);
